import numpy as np
import exceptions

class FractionalPowerPolynomialKernel(object):
    """
    Fractional Power Polynomial Model class (default for PCA ) 
    k(x,y) = sign(xT.y+c)*abs(xT.y+c)^d
    """
    def __init__(self, c=0.0, d=0.8):
        self.c = c
        self.d = d
        if (d >= 1.0) or (d <=0 ):
            raise Exception("d must in range (0.0, 1.0)")

    def __call__(self, x, y):
        dot_result = float(np.dot(x, y)) + self.c
        return np.sign(dot_result) * np.abs(dot_result)**self.d

class Kpca(object) :
    """
    Kernel PCA class 

    implemented very naively and probably unefficiently .. do not use 
    on really big dataset (>few hundreds points) 

    ex : 
    >>> k = Kpca(data) # instanciate 
    >>> k.run()        # compute eigen vectors 
    >>> p = k.project(a_vector)  # get n first projections for any vector 
    >>> p = k.projected_data()   # ....       for the input dataset

    """
    def __init__(self,dataset,kernel=FractionalPowerPolynomialKernel(),components=5):
        """
        @param dataset MxN  : M row samples of dimension N 
        @param kernel is one of the kernel class (must be a callable with two vectors -> scalar ) 
        @param components : # of components keeped for projections 
        """
        self.data = dataset
        self.M = dataset.shape[0]
        self.dim = dataset.shape[1]
        self.kernel = kernel 
        self.comp = components

    def _compute_kmatrix(self):
        """ compute the Kernel dual matrix on the input dataset """
        self.K = np.empty((self.M,self.M))
        for i in range(self.M):
            self.K[i,i] = self.kernel(self.data[i],self.data[i])
            for j in range(i+1,self.M) :
                self.K[i,j] = self.K[j,i] = self.kernel(self.data[i],self.data[j])
                

    def run(self) :
        """
        Compute the eigen vectors of K matrix and store them .. 
        """
        self._compute_kmatrix()
        lamb,a = np.linalg.eigh(self.K) # K is symmetric .. 
        self.alpha = []
        for i in reversed(np.argsort(lamb)[-self.comp:]) : # according to np doc , ev are not necessarily sorted .. 
            self.alpha.append(a[:,i])
        #print self.alpha

    def project(self,x) :
        """
        get the projection of the vector x on the first components 
        """
        res = []
        assert x.shape == (self.dim,)
        for a in self.alpha :
            k = np.array([self.kernel(d,x) for d in self.data])
            res.append(np.dot(a,k))
        return res

    def projected_data(self) :
        """ 
        return the first components of the projected dataset 
        """
        _a = np.array(self.alpha)
        return np.dot(_a,self.K).T

def generate_data(n=100):
    mu = np.array([[1,1],[0,1],[1,0]])
    r = []
    for i in range(3) :
        r.append( np.vstack( [np.random.multivariate_normal(mu[i],np.eye(2)*.02)] for _ in range(n/3)))
    res = np.vstack(r)
    np.random.shuffle(res)
    return res


from matplotlib import pyplot

dat = generate_data(500)

k = Kpca(dat)
k.run()
print k.project(dat[0,:])


delta = 0.05
x = np.arange(np.min(dat[:,0]),np.max(dat[:,0]), delta)
y = np.arange(np.min(dat[:,1]),np.max(dat[:,1]), delta)
X, Y = np.meshgrid(x, y)
print X.shape
print x.shape
print y.shape
Z = [np.empty(X.shape) for _ in range(5)]
for i,y_ in enumerate(y) :
    for j,x_ in enumerate(x) :
        p = k.project(np.array([x_,y_]))
        for c in range(5) :
            Z[c][i,j] = p[c]

print "done projecting"

# Create a simple contour plot with labels using default colors.  The
# inline argument to clabel will control whether the labels are draw
# over the line segments of the contour, removing the lines beneath
# the label

for c in range(5) :
    pyplot.figure()
    pyplot.plot(dat[:,0],dat[:,1],'.')
    CS = pyplot.contour(X, Y, Z[c])
    pyplot.title("comp %d"%c)

pyplot.figure()
p = k.projected_data()
pyplot.plot(p[:,0],p[:,1],'k.') 
pyplot.title("projected data")
pyplot.show()
